Entropy-based measure of structural order in water 
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We analyze the nature of the structural order established in liquid TIP4P water in the framework 
provided by the multi-particle correlation expansion of the statistical entropy. Different regimes are 
mapped onto the phase diagram of the model upon resolving the pair entropy into its translational 
and orientational components. These parameters are used to quantify the relative amounts of 
positional and angular order in a given thermodynamic state, thus allowing a structurally unbiased 
definition of low-density and high-density water. As a result, the structurally anomalous region 
within which both types of order are simultaneously disrupted by an increase of pressure at constant 
temperature is clearly identified through extensive molecular- dynamics simulations. 
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Water undergoes, under compression, a continu- 
ous structural transformation from an open, hydrogen- 
bonded tetrahedral structure to a denser form in which 
the local molecular arrangement beyond the first co- 
ordination shell looks deeply modified 0. Such two 
forms are ordinarily referred to as low-density and high- 
density water (LDW, HDW). There is wide consensus 
on their being high-temperature manifestations of two 
glassy polymorphs of ice: low-density and high-density 
amorphous ice, respectively Q. Raman and Brillouin 
spectroscopic studies have recently given conflicting pre- 
dictions on the thermodynamic boundary between LDW 
and HDW forms 0, 0|. In this Letter we present an 
approach to analyze different ordering regimes in water 
that is entirely based on entropy and which can be used 
to trace the LDW-HDW boundary without resorting to 
any a priori assumptions on the local structure of the 
stable liquid phase. 

Quantifying the degree of structural order present in 
liquids and amorphous materials is a challenging issue 
that has received considerable attention in the last few 
years 0, S 0, 13 • The possibility of characterizing in a 
synthetic way the macroscopic state of a non-crystalline 
substance through some integrated measure of the overall 
amount of order emerging from microscopic correlations 
has been promisingly exploited also for complex molecu- 
lar liquids such as water |9|, [lOl, . The method is essen- 
tially based on the identification of appropriate metrics 
for both positional and angular order. The thermody- 
namic states of the system can then be mapped onto 
a structural diagram, spanned by such two parameters, 
which can be used to obtain a direct and quantitative 
insight into different ordering regimes achieved by ma- 
terials under cooling and/or compression. In this Letter 
we introduce and analyze a metric for discussing struc- 



tural anomalies in water that is entirely and consistently 
based, for both translational and bond-orientational or- 
der, on the "fine structure" of the configurational entropy 
resolved through spatial correlations of increasing order. 

The basic relation between entropy and order can be 
cast in a formally rigorous statistical-mechanical frame- 
work by resorting to the multi-particle correlation expan- 
sion of the configurational entropy 0, 0| . For a classi- 
cal, homogeneous and isotropic fluid, this expansion has 
the form: 



n=2 



(1) 



S2 



{(7(r, in[g(r, w^)] _ g{y.^ _^ IjdrdtJ^ , 



where Sox is the excess entropy per particle in units of 
the Boltzmann constant and the partial entropies s„ are 
obtained from the integrated contributions of spatial cor- 
relations between n-tuples of particles. The pair entropy 
of a molecular fluid can be written as iI4 . J^J : 

2172 
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where p is the number density and g{v,uP') is the pair 
distribution function (PDF) which depends on the rela- 
tive separation r between a pair of molecules and on the 
set of Euler angles = ^ that specify the relative 
orientations of the two particles in the laboratory refer- 
ence frame. The quantity U, represents the integral over 
the Euler angles of one molecule (fi = Stf^ for non-linear 
molecules). The pair entropy typically accounts for the 
overwhelming contribution to the total excess entropy of 
a dense liquid. In order to highlight the relative con- 
tributions of translational and orientational correlations 
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we follow Lazaridis and Karplus [l5j and factorize the 
PDF as 

g{r,uj^)^g{r)g{uj^\r). (4) 

In Eq. 01 g{r) is the radial distribution function (RDF) 
for some arbitrary site in the molecules which, in the 
case of water, is identified with the oxygen atom. This 
function is obtained upon integration of the PDF over 
its angular coordinates. Instead, the function g{uj'^\r) 
represents the conditional distribution function for the 
relative orientation of a pair of molecules whose relative 
distance is r. This quantity is normalized to fl^ and 
is referred to as the orientational distribution function 
(ODF). Using Eq. 0] one obtains for the translational 
and orientational contributions to the pair entropy: 

4*'^--^P l[9{7-)\ngir)-g{r) + l]dr, (5) 

and 

si°'^ ^ p I g{r)S{r)dr , (6) 

with 

Sir) = J g{u;'\r)ln[g{u;^\r)]du;' . (7) 

The use of the positive-definite quantity — ^ mea- 

sure of translational order in both atomic and molecu- 
lar fluids has been already suggested by Truskett and 
coworkers '&] and by Errington and Debenedetti ^9] as 
an alternative to the integral of the absolute value of the 
total correlation function, \g{r) — l|, over a finite number 
of coordination shells. As requested for an order param- 
eter, both quantities vanish in an ideal gas. As for the 
orientational order, the metric exploited so far for water 
is actually a measure of the degree of tetrahedrality in 
the average distribution of the four ox yge n atoms that 
are closest to a given central atom |^ lift fill llfij . As 
such, the definition of this order parameter rests on the 
knowledge of a pre-defined local structure. We propose, 
instead, to adopt —52™^ as an independent measure of 

angular order, au pair with —52*'' for positional order. 
This mutually consistent choice for the translational and 
orientational metrics is rooted on a general and rigorous 
statistical-mechanical ground. As such, it can be im- 
plemented for any molecular fluid, without any a priori 
cognition of the local environment. 

We calculated the translational and orientational pair 

entropies for the TIP4P model of liquid water Ea. Ai- 

though more and more sophisticated classi cal p otentials 
are still being developed d El 113, El El El, the 
TIP4P description is, on the overall, very satisfactory 
when compared with experiments l 2 4| . l25l . and its pre- 
dictive power is well established [23, 123 • 



We carried out molecular-dynamics (MD) calculations, 
implemented through the PINY code |23| on a system of 
108 water molecules in a cubic supercell with periodic 
boundary conditions. The sensitivity of the results to 
the size of the sample was tested with a number of simu- 
lations performed at ambient pressure on a system of 512 
molecules. Long-range forces were computed through a 
particle-mesh Ewald method. The distance cutoff was set 
at 8.5 A. The TIP4P model is based on a rigid-molecule 
description, in that it neglects intramolecular degrees of 
freedom. A time step of 2.5 fs turns out to be sufficient to 
ensure a proper dynamical evolution. Simulation times 
were in the 2.0 ns range. The configurations were stored 
every 1 — 2 ps. The resulting RDFs were calculated over 
1000-2000 different configurations with a spatial resolu- 
tion Ar = 0.05 A. The ODF was calculated at inter- 
vals of 10°. We performed independent simulations in 
the NPT ensemble for pressures comprised between am- 
bient pressure and lOkbar in the 210 — 400 K temper- 
ature range, with a step of 10 K. Lower-density states 
(0.95, 0.90, 0.85, 0.80, 0.75 g/cm^), corresponding to neg- 
ative pressures, were investigated with constant-volume 
simulations carried out over the same temperature range. 
The initial set of positions and velocities of each run was 
taken from the last configuration of a simulation carried 
out at the same pressure, at an immediately higher or 
lower temperature. 

At variance with Sj*'^'' , the calculation of Sj""^-* poses se- 
vere difficulties. Indeed, computing and integrating the 
ODF - which, in principle, depends on up to nine inde- 
pendent variables - directly for a polyatomic molecular 
fluid over an extended range of temperatures and pres- 
sures is a formidable task. In practice, in the case of 
water, the ODF is a function of six variables only: i.e., 
the intermolecular distance, the two angles formed by 
the dipole vector of each molecule with the intermolec- 
ular axis, the angle describing the rotation of the two 
molecules around the intermolecular axis, and the two 
angles describing the rotation of each water molecule 
around its dipole vector. In order to tackle this fairly 
demanding computational task, we resorted to an ap- 
proximate mapping of the ODF onto its low-density 
form. This mapping, known as "adjusted gas-phase" 
(AGP) a ppr oximation, was introduced by Lazaridis and 
Karplus [l5| and is implemented under a two-fold con- 
straint: i) the requirement that the "exact" - i.e., nu- 
merically evaluated - orientationally averaged interaction 
energy between two water molecules be reproduced, as a 
function of their relative distance, in the liquid phase; ii) 
the consistency of the approximate ODF with the one- 
and two-dimensional projections of the exact function 
over the five-dimensional space spanned by the angular 
coordinates. The reliability of the AGP approximation 
has been positively checked through a number of indi- 
cators such as the oxygen-oxygen and oxygen-hydrogen 
site-site distribution functions. An independent direct 
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FIG. 1: (Color online) Ordering map: Translational order 

ftr) 

parameter, — Sj , plotted as a function of the orientational 
order parameter, — Sj""^', for 108 particles at different temper- 
atures (entropy units: cal/mol K). Circles, T = 210 K; stars, 
T = 240 K; squares, T = 270 K; upward triangles, T = 300 K; 
diamonds, T — 320 K; downward triangles, T = 350 K; pen- 
tagons, T — 400 K. Markers are colored according to the 
pressure of the liquid. For T > 300 K, increasing pressures 

(tr) 

lead to increasing values of — Sj , for T < 300 K, increasing 
pressures first produce an increase of both order parameters 
(negative-pressure states), followed by their simultaneous de- 
crease that partially stops at the minimum of the ordering 
locus beyond which the translational order parameter starts 
growing again. The dotted line has vanishing intercept (not 
shown) and was traced as a guide for the eye to follow the 
asymptotic behavior of the order parameters at high temper- 
atures and low pressures. The inset shows a comparison be- 
tween the results obtained at T = 300 K with 108 (triangles) 
and 512 particles (crosses). 



calculation [? ] of the orientational pair entropy of 
512 TIP4P water molecules carried out at constant vol- 
ume (p = 0.999 g/cm^) at T = 298 K leads to a value 
(—52°'^'' = 50.08 Joule/molK, [s^) that is in acceptable 
agreement with the constant-pressure (1 atm) AGP esti- 
mates obtained with 216 molecules at the same temper- 
ature (-5^°''^ = 48.74 Joule/molK, 31]), and with 512 
molecules at T = 300 K (-s^"^ = 48.70 Joule/molK, 
this work). We note that the pair entropy already rep- 
resents more than 80% of the total configurational en- 
tropy of TIP4P water at the highest temperature and 
largest pressure that we have investigated. Moreover, 
this fraction increases with decreasing temperatures and 
pressures until — S2 eventually overcomes — Sex H^- On 
the other hand, orientational correlations account for 
the overwhelming contribution to the pair entropy. We 
found that this contribution also increases, with decreas- 
ing temperatures and pressures, from a percent value 
of - 67% for T = 400 K and P = 10 kbar to - 79%, 
about which it appears to saturate at low enough tem- 
peratures and pressures. We now turn to a discussion of 
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FIG. 2: Phase diagram of the TIP4P model for water: solid 
circles and open triangles identify states corresponding to 
the maxima observed in —s'^'^\p) and —s'^'^\p), respectively; 
open circles identify states corresponding to the minimum ex- 
hibited by —s^'^''{p). Stars locate the states of maximum liq- 
uid density. Lines traced across the data points are a guide 
for the eye. The solid-solid and solid-liquid boundaries are 
sketched after the coexistence lines shown in Fig. 1 of (2^ for 
the TIP4P model. 

the thermodynamic behavior of — Sj''"* and —s'^'K The 
former quantity behaves differently according to temper- 
ature. For T > 300 K, the translational order parameter 
monotonically increases - over the explored range - as a 
function of density or pressure, as is typically observed 
in an ordinary simple fluid Below ambient temper- 
ature, this behavior is still maintained at either low or 
high enough pressures: for T < 295 K, a window actu- 
ally opens up within which — S2*''' decreases with increas- 
ing pressures. As for the orientational order parameter, 
this quantity shows a maximum as a function of den- 
sity, which sharpens with decreasing temperatures. For 
T < 280 K, the maximum observed in both quantities falls 

- to within the numerical uncertainty of the calculations 

- around the same density 0.92 g/cm'^). 

The nature of the correlation between translational 
and orientational order can be enlightened upon con- 
structing a two-parameter ordering map, as originally 
proposed by Torquato and coworkers in The map 
generated by the order parameters introduced in this Let- 
ter is shown in Fig. 1, where we plotted — Sj*''' as a func- 
tion of — Sj"-*. We identify three distinct regimes accord- 
ing to the nature of the structural response of water to 
compression. At high temperatures, both translational 
and orientational order are enhanced by pressure in low- 
density water. However, at higher densities, a different 
structural condition develops, characterized by opposite 
behaviors of the two order parameters: in this hybrid 
regime the translational order keeps increasing while the 
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orientational order decreases as the pressure is raised fur- 
ther and further. At low temperatures, the compression 
of the hquid may also weaken both types of order: such a 
twofold anomalous behavior is observed in an intermedi- 
ate range of densities whose borders on the ordering locus 
coincide with either the relative or absolute maximum of 
— 52°"^'' or — S2*'^\ respectively, and with the minimum of 
— Sj*'''. In passing, we note that the shape of the order- 
ing map is not appreciably modified by the size of the 
calculation (see the inset in Fig. 1). 

The mapping of the (T, P) domains corresponding to 
the different ordering regimes outlined above onto the 
phase diagram of TIP4P water is shown in Fig. 2. The 
lower boundary of the structurally anomalous region is 
entirely located in the negative pressure range. This 
threshold shifts to higher pressures with increasing tem- 
peratures. However, it is interesting to note that this 
line traces, almost all over the range, a constant-density 
path (p ~ 0.92 g/cm^). We surmise that, in the light 
of the present results, this density can be interpreted as 
a sort of intrinsic threshold - in the liquid phase - for 
the nucleation of a stable solid with perfect tetrahedral 
co-ordination. Indeed, the density of the TIP4P Ih ice 
at ambient pressure (0.937 g/cm^) (23i] is consistent with 
this hypothesis. We also recall that the corresponding 
experimental value is ^ 0.92 g/cm'^ [s^j ). 

The upper boundary of the anomalous region shifts 
to lower pressures with increasing temperatures. Corre- 
spondingly, the density decreases from ~ 1.2g/cm'^ to 
^ 0.95 g/cm'^. As seen from the figure, a twofold anoma- 
lous structural behavior is observed in the liquid coexist- 
ing with ices Ih, III, and, partially, V. It is rather natu- 
ral to explain such an anomalous behavior as due to the 
presence of a persistent and extended hydrogen-bonded 
network. Consequently, we are lead to identify the asso- 
ciated region as the LDW structural basin, as opposed to 
the HDW basin where the translational and orientational 
order parameters show opposite trends as a function of 
the pressure. The lower boundary of the HDW region, 
separating above ambient temperature thermodynamic 
regions where the angular order is either enhanced or 
disrupted by pressure, crosses the P = axis at a tem- 
perature (r ~ 350 K) that is close to the boiling point 
of TIP4P water at ambient pressure (T ~ 363 K) [33 |. 
We also note that the locus of density maxima lies en- 
tirely inside this region as also found by Errington and 
Debenedetti through a tetrahedricity-based order repre- 
sentation 0. 

On the basis of the analysis of the experimental data 
carried out by Soper and Ricci 0, one is lead to iden- 
tify the LDW-HDW crossover threshold with the state 
of the liquid whose oxygen-oxygen structure factor can 
be represented by an equal-weight superposition of the 
"asymptotic" LDW and HDW forms. At T = 268 K, 
this criterion yields a pressure of 1 kbar. The esti- 



mate obtained with the present method for the TIP4P 
model at the same temperature is about 1.2 kbar. As for 
the thermodynamic slope of the LDW-HDW structural 
boundary, our results are consistent with the numerical 
predictions of Saitta and Datchi who used a cri- 
terion based on bond-angles distributions to trace the 
crossover line, and with the experimental findings of Li 
and coworkers 0. 

In this Letter we have implemented, on a rigorous 
statistical-mechanical basis, an entropy-based theoretical 
framework that allows a structurally unbiased descrip- 
tion of ordering regimes in liquid water, in its thermo- 
dynamically stable domain. The current individuation 
of LDW and HDW basins is in good agreement with 
neutron-diffraction results at 268 K and further clarifies 
other conflicting experimental data. We plan to extend 
this analysis to the study of amorphous ices. However, 
the proposed metrics prospectively qualify the current 
method as a general tool for investigating the nature of 
structural order in molecular liquids. 
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